Motor effect in electron transport through a molecular junction with torsional 

vibrations 
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We propose a model for a molecular junction with internal anharmonic torsional vibrations in- 
teracting with an electric current. The Wangsness-Bloch-Redfield master equation approach is used 
to determine the stationary reduced density matrix of the molecule. The dependence of the cur- 
rent, excitation energy and angular momentum of the junction on the applied voltage is studied. 
Negative differential conductance is observed in the current-voltage characteristics, ft is shown that 
a model with vibrationally dependent coupling to the electrodes, asymmetric with respect to the 
interchanging of electrodes, leads to a strong correlation between the applied voltage and the an- 
gular momentum of the junction. The model thus works as a molecular motor, with the angular 
momentum controlled by the size and sign of the voltage. 

PACS numbers: 73.23.-b, 71.38.-k, 85.65.+h, 85.85.-l-j 
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I. INTRODUCTION 

Molecular electronics has become a dynamically grow- 
ing field in the last decade. It is a highly interdisciplinary 
research topic presenting many challenges in both the- 
oretical and synthetic chemistry, solid state and many 
particle physics and non-equilibrium statistical physics. 
Much attention has been paid to the careful prepara- 
tion of single molecule junctions and to studies of the 
conductive properties of different molecules. Numerous 
papersir— , studying both the experimental and theoret- 
ical aspects, have already been devoted to the subject. 
It has also been notedSr— that the molecular structure 
and vibrations or conformational changes are among the 
main points of interest, making the molecules distinct 
from other electronic elements. 

The coupling of electronic and mechanical degrees of 
freedom is a standard part of electrical engineering and 
the resulting gadgets are a standard part of our life. In 
recent decades the coupling of electronic and vibrational 
degrees of freedom has been achieved in nanoscale solid 
state devices (see Schwab and Roukes^ and Blencowe-^ 
for reviews on NEMS) with the quantum regime reached 
in the mechanical degree of freedom. Driving the molecu- 
lar vibrations with the electronic current is a natural ex- 
tension of this concept. We can consider a molecular elec- 
tronic element or molecular junction with some molecu- 
lar groups performing rotational motion in response to 
the bias voltage across the molecular junction. Such el- 
ements have already been anticipatedii and the work is 
closely related to electron shuttle o^^'^'^ . Recently, the 
excitation of periodic nuclear motion in molecular junc- 
tions has been studied as a classical motion of atoms in 
non-conservative forces induced by the current flow^^. 
Artificially-built molecular motors, where molecular vi- 
brations are driven by light or stochastic fluctuations 
due to interaction with a thermal bath (Brownian mo- 
tion), have also been studied. The externally-driven tor- 



sional motion of some small parts of molecules has been 
demonstrated for molecules both in gas-- and mounted 
on surfaces^. 

The main goal of this paper is to demonstrate that 
the rotational motion of a molecular group in a metal- 
molecule-metal junction should in fact be a very common 
phenomenon and that there are only two conditions re- 
quired: 1) the presence of some part of the molecule capa- 
ble of rotation with a moderately small potential barrier 
against this rotation; and 2) a breaking of the mirror 
symmetry in the junction (chirality of the junction). 

To achieve this goal we set up a general model for the 
description of the interaction of a current flowing through 
a molecule with anharmonic molecular vibrations. To an- 
alyze and understand the effect of anharmonicity in more- 
or-less well-controlled conditions we first define a model 
with some small vibrational coupling due to a small vi- 
brational potential energy shift. We then switch to a 
more realistic model of molecular vibrations motivated 
by real molecular rotors as used in a previous experiment 
with light-driven artificial molecular motors^^. The dy- 
namics of the system will be studied using rate equations 
for the reduced density matrix of the molecule. The cur- 
rent through the junction and the average angular mo- 
mentum of the molecule are calculated as a function of 
the voltage drop across the junction. 



II. MODEL 

The molecular junction consists of two metallic elec- 
trodes or leads L and a molecular bridge M connected 
between them. The corresponding division of the Hamil- 
tonian reads 



H = Hm + H, 



H 



ML- 



(1) 



The model of the bridging molecule is assumed to con- 
sist of one vacant electronic level that allows electrons 
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FIG. 1. Schematic representation of the junction 



to tunnel through it from the left lead to the right lead. 
We will call this state the localized state. There is an 
extra charge on the molecule when this state is occupied, 
and we will thus speak about a neutral, or a charged, 
molecule (anion) if the state is unoccupied or occupied, 
respectively. In addition, we include a vibrational degree 
of freedom on the bridge, which exchange energy with 
the electrons. We thus assume 



Hm ~ h^dd} + hi<v d, 



(2) 



where d' and d create and annihilate electrons in the lo- 
calized state on the bridge. The operators h^ and hi 
describe the vibrational motion of the nuclei in the neu- 
tral and charged state respectively. The form of both 
ho and hi is based on the Born-Oppenheimer vibra- 
tions for isolated molecules and in the spirit of the Born- 
Oppenheimer approximation we assume, that both /iq 
and hi commute with d^ and d. The Born-Oppenheimer 
approximation is, of course, broken when we allow elec- 
tronic transitions due to the coupling to the leads Hml- 
Figure [U showing schematically the system that we 
have in mind, may serve to guide us through further 
specifications of the model. We will consider only one 
nuclear (vibrational) degree of freedom representing the 
rotation of a part of the molecule, pictured as a ben- 
zene ring bound in para positions to other parts of the 
molecule that are directly attached to the metal leads. 
The vibrational Hamiltonians /iq and hi are both of the 
form 



' 21 9(^2 



^.(^), 



(3) 



where tp £ (0, 27r) represents the vibrational coordinate 
(angle of rotation of the ring) , / is the moment of inertia 
of the ring and Vi{ip) is the Born-Oppenheimer vibra- 
tional potential of the molecule with unoccupied (i = 0) 
or occupied (i = 1) electronic level. In this paper we 
want to characterize the main features of the dynam- 
ics of transport through junctions with such anharmonic 
vibrations and we assume a simple analytic yet rather 
general shape of potentials'^ as 



Vi{(p) = e.j + Ai cos(nj</3 -t- ipi). 



(4) 



We use two sets of models. In the models of group 1 we 
want to capture the basic features of the usual models 



used for studing molecular conduction junctions coupled 
to vibrationsi?:"— , but to allow for large amplitude an- 
harmonic motion. We thus take both Vq and Vi to be the 
potentials of the mathematical pendulum (no — ni — 1), 
one shifted with respect to the other (eq, ipo differs from 
El, ipi). For these models we also set the moment of iner- 
tia / to an unrealistically small value-- to allow for a more 
efficient numerical solution. The second group of mod- 
els are motivated by more realistic parameters, as might 
be expected in real molecular systems (see, for example, 
Fortrie and Chermette— or Tsuzuki et al.—). The vi- 
brational potential Vq for the neutral molecule is thus 
characterized with a smaller amplitude Aq and larger 
number of oscillations (uq — 2) than the potential for 
the charged molecule Vi (for biphenyl junction see Ci'zek 
et al."^^). Also the value of / is set larger to describe 
the inertia of the benzene ring. All values of the model 
parameters are summarized in Table ID 

In principle, the parameters of the molecular Hamilto- 
nian should also depend on the voltage applied across the 
junction. We follow the common practice (see, for exam- 
ple, Galperin et al.— or Hartle et al.— and the works 
cited there) of using the molecular Hamiltonian indepen- 
dent of voltage. This approximation is also supported by 
the calculation of molecular potentials in the electric field 
made by Petreska et al.^-, where the potential barriers 
do not show a significant change for realistic values of the 
electric field. 

The Hamiltonian of the leads is written in the form 



Hr 



7 . 7 ,£kaC^.^Cka, 
a—l.r k 



(5) 



where the operator cj,^ creates an electron in the state 
with a wave number k in the lead a S {l-,r}. To be 
more specific we assume that the Hamiltonian ([5]) was 
obtained from the diagonalization of the one-dimensional 
nearest-neighbor tight-binding model^i, i. e. the disper- 
sion relation reads 



Ska = A*a +2/3iCOs(fc), 



(6) 



where fia is the chemical potential of the lead a. When 
we talk about the voltage U applied to the junction we 
assume m = +^ for the left lead and fir = —^ for 
the right lead. The parameter (3i = 3eV defines the 
width of the conduction band in the leads. It is chosen 
in order to provide a band, wide enough to exclude edge 
effects. Each lead is separately assumed to be in thermo- 
dynamic equilibrium with states populated according to 
the Fermi-Dirac distribution 



fai^ka) — 1 



=,(efec-Mc)/fcT 



(7) 



To complete the description of the model the molecule- 
lead coupling 

Hml = ^Y. ^dkai^^'cka + cl^d) (8) 

a—l.r k 



TABLE I. Summary of all parameters for the models in use. The energies Si, Ai are in units of eV, angles in radians and 
moment of inertia in atomic units (m^aQ). 
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must be specified. We assume the separable form, of 
the coupling coefficients Vdka — VkVaiv)- The depen- 
dence on the electron momentum Vk — (^2 sin(fc) is again 
motivated by the one dimensional nearest-neighbor tight 
binding model, where the sine term comes from the elec- 
tronic wavefunction^S and the parameter P2 is the overall 
coupling strength, which is set to P2 = 0.07 eV. We shall 
study several forms of the angle dependent part Va((p). 

• Model la. The case of angle- independent coupling 
Va{^p) = 1. This assumption simplifies the treat- 
ment of the dynamics significantly. Furthermore, 
this is the case most often considered by other 
studieai^".2i. 

• Model lb. Here we introduce angle- dependence into 
the coupling term Va{^p) = cos((/3 — Lpa) but we 
do not break the symmetry between the left and 
right lead tpi — (fir = tt. The form of the angular 
dependence can be motivated by the Hiickel model 
(see Cizek et al.'^^ , Pauly et al.^^). 

• Model Ic. Finally we break the symmetry of the 
system, taking Va {^p) in the same form as in model 
lb but with Lpi = iT^ ipr = 37r/2. The asymmetry 
can arise as a result of different molecular bonding 
to the left and right lead. We can, for example, 
consider a molecule consisting of a chain of three 
aromatic rings with the first ring fixed to the left 
lead the middle ring acting as the rotor, and the 
right ring fixed to the right lead (this idea is used 
in the diagram of the model in Figure [Ij . Breaking 
the symmetry provides circumstances for observing 
the " motor effect" , i.e. preferential rotation of the 
rotor in one direction depending on the sign of the 
voltage applied across the junction. 

• Models 2a and 2b use the same form of coupling as 
the Model Ic. The parameters ipi, ipr for all models 
are summarized in Table ID 



set of projectors 



dd^ +dU^l, 



(9) 



where the projector dd^ projects on the part of the 
Hilbert space with the unoccupied molecular bridge and 
d^d projects on the occupied bridge space. It is advan- 
tageous to use different basis sets in the vibrational part 
of the Hilbert space according to the occupation of the 
bridge. We thus define two basis sets \n) and \v) 



ho\n} = En\n), 
hi\v) — E^\v). 



(10) 



This choice diagonalizes the molecular part of the Liou- 
villian operator Cq (see next section) and in the case of 
the harmonic molecular potentials Vi it is equivalent to 
performing the polaron transform^S as made in the in- 
dependent boson modeP". With cosinusoidal potentials, 
the states |n), \v) can be expressed in terms of Math- 
ieu functions. Energy levels En and E^ (n,v = 0,1,...) 
of the unoccupied and occupied molecule for models 1 
and 2 are shown in Figures [5] and [3] respectively. From 
the point of view of classical mechanics, such a shape of 
the potential provides a rotational barrier with an en- 
ergy of ei -\- Ai, which divides the two types of motion: 
vibrational (when the energy of the system is below the 
barrier) and rotational. In a quantum mechanical de- 
scription the wavefunctions of states with an energy be- 
low the barrier are localized in space and in this sense 
can be called vibrational. In contrast, states above the 
barrier are delocalized through the whole interval (0, 2n). 
These states are two times degenerate, as the two direc- 
tions of rotation are available above the barrier and they 
are energetically equivalent. We will call them rotational 
states. Well above the classical barrier, when the sys- 
tem has a lot of energy and does not feel the potential 
anymore, the states almost coincide with the free rotor 
states. 



III. THEORY 

The Hilbert space of our system is a direct product of 
the spaces of electronic and vibrational degrees of free- 
dom. In the electronic space, we can define a complete 



A. Master Equation 

Different theoretical approaches to the transport of 
charge across a molecular junction with a coupling to 



where the Liouvilhan super operator £ = £o + £i with 




Angle [rad] 



FIG. 2. Energy levels of the unoccupied (red lines) and occu- 
pied (green lines) molecular bridge of Model 1. The respective 
potentials are the borders of dark gray and light gray areas. 




Angle [rad] 

FIG. 3. The energy levels and potentials for Model 2 plotted 
in the same way as in the previous figure. The angular de- 
pendence of molecule- lead couplings Vaif) = cos(iy3 — ipa) are 
shown in the inset for the left a = I and right a = r leads. 



Ci[p{t)] = -Tl-L JdT[HML, [HML{~T),p{t)(g)pl]], 


(12) 
and 

Hml{-t) = e-'(^«+^^)"i7MLe+'(^«+^^)". (13) 
These equations are derived as a second order expansion 
in Hml, which is assumed to be small. The equilibrium 
RDM of the leads is denoted p^. The first term Cq de- 
fines the time evolution for a " free" molecular bridge that 
is disconnected from the leads. The second term Ci in- 
corporates the influence of the leads and could be written 
as the sum of two independent terms £i — Cij + ^i,r for 
the left and right leads respectively. 

Here we are not interested in the time evolution of the 
density matrix. We assume the existence of a stationary 
state dp/dt = which satisfies the equation C[p] =0. In 
basis representation, the Liouvillian is a 4th rank tensor. 
We search for a nontrivial solution of the equation 



2_^^ijklPkl — 0. 



(14) 



fc; 



Before we write expressions for the components of this 
tensor we decompose RDM p in the following way 



p = poodd^ + piid^d + poid + p^d^ . 



(15) 



With this representation the right hand side of the equa- 
tion pT|) can be reorganized in blocks 



(16) 




vibrations have been discussed in the review article of 
Galperin et al.—. Here we consider a weak-coupling case, 
which can be treated with the standard master-equation 
approacb^^^i^. A slight modification is needed to ac- 
count for the anharmonicity of the vibrations as de- 
scribed in what follows. 

In the framework of master equation (ME) theory 
we calculate the reduced density matrix (RDM) of the 
molecular bridge p. We start from ME in the Wangsness- 
Bloch-Redfield (WBR) formes 



Equations for poi and pio are decoupled from this system 
for Poo and pn and they are not included here, since the 
observables of interest are also independent of poi and pio 
(see below). To write the elements explicitly we use basis 
sets pO]) . Basis {\n)} for the block p„„/ = (n|poo|'^') and 
{\v)} for the block p^j^i = (u|pii|i;'). We will consistently 
use the letters n, rii, 77.2, n' to number the vibrational 
states of the unoccupied molecule and the letters t;, t;i, 
V2^ v' for the states of the occupied molecule, omitting 
the index 0/1 that would otherwise distinguish the occu- 
pancy. Thus the components /:°-„-„^„^, ^y'^^'^^,^, of the 
zeroth order contributions Cqq, /Z^i read 



V[,V'^,VI,V2 






(17) 



-p{t) = cMt)] 



A[p(<)], 



The lowest order contributions £"■ describing the pres- 
(11) ence of the leads a = Z, r in eq. p6|) are expressed as 



n-^ ,n2,ni .122 

fOi 

n'^ ,n2,vi ,f 2 

fOi 



V "^ V ^ 

■n. ^ n ^ 



(18) 



where a;„ 



Ey — En is the transition energy, V"^ 



{n\Va{ip)\v) is the generahzed Franck-Condon overlap 
for the transition and the factors fa{(^)Ta{tjj) and [1 — 
fa(LL!)]Ta{uj) come from the imaginary part of the func- 
tions 



UE) = iJdTj:e 



-i(eka-E)T 



ifkcWl 



UE) ^ I J dr E e+^(^--^)^(l - fkcW^ 



(19) 



resulting from the time integration in eq. (J12p . The real 
part of these functions that leads to renormalization of 
the energy levels is neglected here (see also discussion in 
Leijnse and Wegewijs'^^ where it is argued that the real 
part is canceled out in higher orders). The imaginary 
parts can be calculated analytically. For the tight binding 
model of leads T{E) become^*' 



/3| 



r{E) = ^^^4(3i^iE-f,^y 



(20) 



inside the band (i.e. when E € [/i — 2/3i,/i + 2/3i]) and 
equal to zero outside. 

It is often argued that the nondiagonal elements (pnn' 
and pyyi for n ^ n', v ^ v' ) in RDM decay rapidly in 
time and are neglected in the search for the stationary 
state. Here the nondiagonal elements for near-degenerate 
states lead to a nonzero angular momentum for the ring 
and can not be neglected. However, to increase the nu- 
merical efficiency we consider only the elements close to 
the diagonal and assume that the RDM has a band struc- 
ture. Final results are presented for the RDM, which 
includes 14 sub diagonals, while its overall dimensions 
are 202x202 for Model 1 and 402x402 for Model 2. We 
tested that the presented results are stable with respect 
to changes in the number of sub diagonals and the num- 
ber of basis functions used in the calculation. 

For numerical convenience, we reshape the tensors C 
of the fourth rank into matrices using compound indices, 
mapping the pair of numbers nn' to a single number v 
(and similarly for vv'). Instead of p4)) we then solve the 
set of equations 



'^C^„'P„' = 



(21) 



together with the normalization condition Tr{p} = 
J2n Pnn + J2v Pw = 1 to find the stationary state. 



B. Observables of interest 

The general formula for the current (see, for example, 
Hartle et al. ^^ ) through one level in ME theory reads 

/ = / drTr i [Hml{-t), P ® P?] ^ Vdk{d^ Ck - eld) \ , 

Q K k ) 

(22) 
where only the contribution from the left lead is included 
in Hml- In our case, the formula reduces to 



nv \ 



'l{uJnv)y^Vj^,yPn'n' 



il-fl{uJnv)]J2Kv'P^^' \- 
v' ) 



(23) 



The mean value of the molecular Hamiltonian Hm gives 
the average excitation energy of the bridge. This quan- 
tity can be expressed as the sum of two contributions 
from the molecule occupied/unoccupied with an addi- 
tional electron, and in similar fashion as was the current 

{Hm) = (£o) + {El) = V Enpnn + V E,p,y. (24) 



The last important observable to be discussed is the an- 
gular momentum of the molecule. This is by construction 
constrained along the rotational molecular axis z in our 
model. The corresponding operator reads 



OLp 



(25) 



We notice, that the operator L^ acts independently on 
the occupied and unoccupied bridge spaces, or, in other 
words, the off-diagonal blocks poi ^-nd pio in the expan- 
sion (ITS]) do not contribute to the mean value. It allows 
us to write the mean value in the form 



(L^) = ^{n\L^\n')pn,n + '^{v\L^\v')py'y 



= 2^^Pn 






mpl^, (26) 



where we have defined the populations p^, p^ of ro- 
tational eigenstates Lz\m) = m\m) for unoccupied and 
occupied states of the molecule respectively, i. e. 



Pm = ^(m|n')p„'„(n|m), 

7171' 

pin ^y^(m\'"')pv'v{v\m). 



(27) 
(28) 



Notice the importance of the off-diagonal elements of 
the RDM in (pS)) . Since the diagonal elements {n\Lz\n), 
{v\Lz\v) are zero we get (L^) = if we neglect the off- 
diagonal elements p„„', Pw' for n y^ n', v ^ v'. While 
the value of the current is mainly determined by the di- 
agonal elements /?„„, pw the nonzero angular momen- 
tum is a consequence of the non-vanishing off-diagonal 
elements in the energy representation, or equivalently a 
consequence of the asymmetry pl^ ^ pL^ in the angular 
momentum representation. 



IV. RESULTS AND DISCUSSION 

In this section we discuss the results for the calculation 
of the current- voltage characteristics and other proper- 
ties of the junctions. We start with Model 1, where the 
harmonic approximation holds for the small amplitude 
of vibrational motion and the levels are near equidistant. 
In addition to the known behavior of such junctions with 
harmonic vibrations we can expect some new effects due 
to the dependence of the molecule-lead coupling on the 
vibrational coordinate. For any higher vibrational excita- 
tion of the junction molecule we would expect a breaking 
of the harmonic approximation. 



A. Current and excitation function 

The current- voltage characteristics and the vibrational 
excitation energy of the junction for Models la-c at tem- 
perature T = 50K are shown in Figure ID Let us first fo- 
cus on the current-voltage curve. The red curve, for the 
model la, exhibits the behavior expected for the model 
with a very small coupling between the vibrations and 
electronic motion. We observe a resonance step at the 
voltage of 0.2 V corresponding to twice the charging en- 
ergy of the molecule 0.1 eV. 

The models lb (green line) and Ic (blue line) differ by 
having a more complicated step structure and the dif- 
ferent stationary value of current that is finally reached. 
This second difference is easily understood: while the 
coupling to the leads has the same maximum strength 
for all of the models (they all reach the value 0.07 eV), 
the angular dependence in models lb, Ic makes it effec- 
tively smaller (a quantitative argument as to why this 
difference is given by factor of two follows at the end of 
this section). The angular dependence is also a source of 
the negative differential conductance since the coupling 
in Model lb reaches a maximum for angles near the equi- 
librium position of the vibrational coordinate ip. The 
non-equilibrium vibrational distribution for larger volt- 
ages therefore reduces the coupling. The negative differ- 
ential conductance effect disappears in Model Ic because 
the coupling to the right lead Vii{(p) peaks at a different 
angle and is therefore effectively increased by the increase 
in the vibrational excitation of the molecule. Both Mod- 
els lb and Ic have the same average value of couplings 




-0.5 0.0 0.5 

Voltage [V] 

FIG. 4. Current (top) and mean value of excitation energy 
(Hm) (bottom) as a function of the voltage, applied to the 
junction, for models la-c at temperature T = 50 K. The 
shaded bars show the positions of the steps derived from the 
energies of the molecular levels (see text). 



yj,/fl((/?) and the asymptotic value of the current for large 
voltages is therefore identical for both models. 

To explain the details of the step-like behavior of the 
curves we start with the mechanism of sequential tun- 
neling through the bridge. Thus electron conduction is 
understood as a sequence of charging (electron attach- 
ment) 



e" -I-M(n) -^ M'{v) 
and discharging (electron detachment) 

M^{v) -> M{n) + e" 



(29) 



(30) 



events on the bridge, where e~ is the electron in the leads 
and M{n) and M~{v) stands for the neutral molecule 
and the anion with the vibrational states \n) and \v) 
respectively. In the first event (charging), the electron 
starts in one of the leads in state |fc) and jumps into the 
unoccupied bridge in the vibrational state \n) and turns 
it into an occupied bridge in the vibrational state \v). 
In the second event the electron starts in the occupied 
bridge in the state \v) and leaves the bridge in the state 
\n) jumping into the leads to the state \k). The energy 



TABLE II. Inelastic one-electron attachment and detachment processes. The table gives the threshold electron energies of the 
inelastic processes and approximate size of the matrix elements responsible for the transitions due to electrons from the left 
and right leads. All values are in units of eV. 



Series 



Process 



^Nv 



Ev — En 



\{n\Vi\v) 



\{n\Vr\v)\' 



El 



ek + KC O Ml^n 



efc + M° « Ml^n+i 
efe + M° ^ Ml^n+2 



0.1 



0.2- 1 



10" 



Exl 
Ex2 
Ex3 



ek + M°^ Ml^n-2 



0.21-0.23 
0.31-0.36 

0.42-0.48 



0.005 - 0.015 

0.001 - 0.06 

< 0.005 



0.05-0.3 

0.001 - 0.01 

< 0.005 



Dxl 
Dx2 
Dx3 



-(0.01-0.03) 
-(0.16-0.12) 
-(0.23-0.28) 



0.005-0.015 

0.001 - 0.05 

< 0.005 



0.05 - 0.3 

0.001 - 0.01 

< 0.005 



is conserved in both of the events 



£k + En 



En, 



(31) 



where Sk is the energy of an electron in a state \k). The 
vibrational energies En and E^ of the occupied and un- 
occupied molecular bridge were defined in equation ([TOl) . 
Each difference uJnv = Ey— En defines a threshold energy 
for one possible in/out channel, which becomes available 
at the voltage U — ±2a;„t, (when the chemical potentials 
of the leads are equal /Zq, = ±-j = ±a;„„) and can (but 
not necessary will) show itself as a step in Figure S) Since 
the structure of the energy levels is the same for all the 
models la-c, we can expect the steps at the same volt- 
age in all three current-voltage curves. The values of the 
threshold energies ujnv = E^ — En are shown in Table HH 
The transitions are divided into several groups. The 
vibrational state of the molecule is unchanged and v = n 
in the first group denoted El. The corresponding thresh- 
old energy is w„,y = O.lel/ independent of the value of n 
because the shape of the two potentials Vo{(p) and Vi{(p) 
is identical except for a vertical and horizontal shift. The 
excitation groups Exl-Ex3 correspond to u = n + 1, 

V = n + 2 and v = n + 3. In the harmonic approximation 
the threshold energies would also be independent of n. 
Table HI] shows the range of values of ujnv for low states 
n = 0, 1, ..., 10 (the choice of this maximum value of n 
is guided by the average excitation energy shown in the 
lower part of the figure) . Lastly, the deexcitation groups 
Dxl-Dx3 are characterized hyv — n — 1, v = n — 2 and 

V = n — 3. If we compare the predicted positions of the 
steps U = zt2ujnv given by the values from Table |TT] with 
the positions of the steps in Figure HI we see that the 
steps correspond to voltages approximately 0.2, 0.3, 0.45 
and 0.7 volts, i. e. to the series El, Dx2, Exl and Ex2. 
It is remarkable that the vibrational energy of the bridge 
doesn't show any steps for a voltage \U\ > 0.8 V and 
grows parabolically for all models. The current-voltage 
characteristics also become smooth in this voltage region. 
The level of excitation of the molecule is too high for the 
harmonic approximation to apply and many different vi- 
brational states are involved. The energy differences be- 
come irregular and it is not possible to sort them into 
groups, as in the case of the low, near harmonic vibra- 




0.04 
0.02 





FIG. 5. Angle distributions and populations calculated inde- 
pendently for the unoccupied and occupied bridge in Model 
Ic. 



tions above. We would expect that quasi-classical theory 
can be applied in this regime. 

At equilibrium {U = 0) only one lowest vibrational 
level En = 0.065 eV is populated (it is localized close 
to the minimum in the vibrational excitation curve in 
Fig. 01). The bridge remains "closed" for current until 
the first step occurs at the voltage 0.2 V, at which point 
the whole El set of channels then opens. If we consider 
the bridge originally in its ground state with energy _EjJ 
the first tunneling event 



£k+En 



El 



El 



>K, 



v=0 Dxl 



> En=0 + Ek 
-> En=i + Ek' 



(32) 



can leave the bridge in the excited state En=i through 
the channel Dxl. When the next electron comes and the 
bridge is already in an excited state, tunneling can excite 
it even higher, but only by one quantum per tunneling 
event. Excitation to higher levels is limited by the com- 
peting process of de-excitation. This picture is the same 
for all three models. To appreciate the differences among 
the models we must look at the Franck-Condon factors 
(n|V;|u)p and |(n|yr|i')P responsible for the strength of 
each ek + M^ -H* M^ transition. 
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FIG. 6. Angular probability distribution of the bridge p{'p) 
at different voltages. 



It is V; = Vr =constant for the Model la. Since the 
potentials Vq and Vi differ only slightly, it is {n\v) ~ 
6nv and the El channels are dominant with only a small 
contribution from Exl, Dxl and virtually no contribution 
for higher channels. This explains vi^hy the red curves in 
the current-voltage and excitation graphs show steps only 
at voltages corresponding to these channels. 

The values of |(n|VQ|w)p for Models lb and Ic are 
shown in Table HIl We should keep in mind that Vi=Vr 
for Model lb. The values of KnlT^I"!^)!^ shown in the last 
column are the ones for Model Ic. There is a pronounced 
difference between Models lb and Ic in the first step at 
U = 0.2 V and we can now understand why. While both 
charging and discharging of the bridge proceeds domi- 
nantly through the El channel for Model lb (giving a 
step similar to that in Model la), the discharging to the 
right electrode through the channel El is strongly sup- 
pressed for Model Ic. Thus for Model Ic discharging of 
the molecule to the right lead is possible mainly through 
the channel Dxl. This gives a smaller value for the cur- 
rent, but a higher value of vibrational excitation for this 
model. Also, for higher voltages, the sizes of the steps 
in the excitation curves follow the sizes of the Franck- 
Condon factors. 

One striking feature apparent from Figure 0] is the 
asymmetry of the curves for Model Ic - a consequence 
of the asymmetry |(n|V;|w)p ^ |(n|K-|w)p. For negative 
voltages the channel El is not available in Model Ic, even 
for [/ < —0.2 V, because charging has to proceed from the 
right electrode and the Franck-Condon factor |(n|K-|ti)p 
is suppressed by four orders of magnitude. Charging of 
the bridge only becomes possible with the availability of 
the Exl process. Both current and vibrational excitation 
is thus only significant for negative voltages U < —0.45 V. 
Another way to look at this asymmetry is to inspect the 
angle distributions 

Poi^)^TT{dd^p\if){^\} (33) 

p^{^) = TY{dUp\^){^\} (34) 

and the populations pyy of vibrational levels on the oc- 



cupied and pnn on the unoccupied bridge, respectively. 
These are shown in Figure [5] for voltages U = ±0.5 V. 
We can make the observation that the angular distribu- 
tion for the occupied molecule follows the shape of the 
angular dependence of the coupling to the donor elec- 
trode (left/right for positive/negative voltage), and the 
distribution of angles for the unoccupied molecule follows 
the coupling to the acceptor electrode. The same effect 
is responsible for the asymmetry found in the population 
distributions in the lower part of Figure [SJ Similar dis- 
tributions for Models la and lb (not shown here) exhibit 
no asymmetry for the change U -^ —U. For these two 
models the vibrational distribution is also weakly corre- 
lated with the charging state of the molecule, i. e. the 
distributions p„„ and pyy have almost identical shape as 
functions of energy En and Ey (i. e. the red and green 
curves in Figure [5] are almost overlapping for Models la 
and lb). The small difference in the populations is only 
due to the small difference between the potentials Vq and 
Vi. The role of symmetries is further discussed in the 
next section. 

The angular distributions p{ip) — TT{p\ip)(ip\} at the 
voltages U — 0.25, 0.5, 0.75 and 1.0 volts are compared in 
Figure[S]for all the models la-c. There is a common trend 
for the angle to become more and more delocalized while 
the voltage grows. The degree of delocalization follows 
from the excitation curve in Figure [5] (lower graph) . At 
U = 0.25 V, the distribution for Model Ic is the broad- 
est, while Model la wins out at higher voltages. The 
angle is completely delocalized above the last step of the 
excitation in the IV curves. The current for Models lb 
and Ic is asymptotically a factor of two smaller than the 
current for Model la. We have just seen that the angle 
distribution at the bridge is more or less homogeneous 
for large voltages. For this reason, the angle-dependent 
part of the coupling for Models lb and Ic, which is equal 
to {cos{(p — fa))'^, reaches its mean value 0.5. For Model 
la, the angle-dependent part of the coupling is constant 
and equal to 1, i.e. twice as large. 



B. Angular momentum (motor effect) 

In Figure [71 we plot the mean value of the angular 
momentum (L^) against the voltage applied across the 
junction. The calculated value of the angular momentum 
(Lz) is in general nonzero for all the models la-lc and can 
reach values of the order of the reduced Planck constant H 
(atomic unit of angular momentum) , the value depending 
strongly on the voltage. The model is thus an example 
of a molecular scale device which can perform rotations 
controlled by the voltage applied to the device. This 
effect is particularly pronounced in Model Ic, where the 
direction of rotation is inversed with the inversion of sign 
of the voltage. 

The symmetries of the curves can be understood from 
the symmetries of each model. None of the models ex- 
hibit symmetry with respect to the inversion of the an- 
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FIG. 7. Mean value of the angular momentum of the molecule 
{Lz) as a function of voltage for Models la, lb and Ic at the 
temperature T = 50 K 
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FIG. 8. Populations p^ and p„ in the angular momentum 
basis for Model Ic at different voltages. 



gle ip -T- — (p, because this change inverts the horizontal 
shift of potential Vi{if) relative to Vb(v3). Breaking of 
this symmetry allows for nonzero values of angular mo- 
mentum in all the models. But Models la and lb are 
symmetric with respect to the mutual interchange of the 
left and right lead. This is reflected in the symmetry of 
the angular momentum with respect to the voltage inver- 
sion U -^ —U. Breaking of this symmetry makes Model 
Ic distinct. 

In all the models la-lc, the molecule does not rotate 
when the absolute value of the voltage is smaller than 
0.4 V. The onset of each curve follows a degree of vibra- 
tional excitation of the molecule (examine lower graph 
of Figure S]). This behavior reflects the fact that the 
molecule should be excited sufficiently high and overcome 
the rotational barrier at 2.5-2.6 eV in order to perform 
rotational movement. A significant population of states 
above the barrier can only be expected when the mean 
value of the vibrational energy of the bridge molecule 
(lower graph, Figure U) is of the order of magnitude of 
leV. The role of these "over-the-barrier" states was also 
checked by omitting all of the states below the barrier 
from the calculation. The linear growth of the angu- 
lar momentum in Model Ic is hardly sensitive to this 
change. Another way to restate this discussion is to look 
again at the angle distribution functions in Figured The 
rotational motion of the molecule is indicated in the de- 
localization of the angle distribution through the whole 
interval (0,27r). At 0.5V, only Model la has this prop- 
erty, but at 0.75 V all three models can rotate. These 
observations are in correspondence with Figure [T) 

The off-diagonal elements of the RDM are often ne- 
glected in the master equation approach, resulting in a 
set of equations only for the population /C»„„, pw of the 
states. In fact, we used this line of thinking when dis- 
cussing the current-voltage and excitation curves in the 
previous section. If we want to capture the motor effect 
we have to consider the off-diagonal elements - at least 
for the near degenerate states close to or above the rota- 



tional barrier. This matter we discussed already in the 
Theory section below Formula ((26)) . where we also saw 
that the coefficients p^ and p^ give a better insight into 
the calculation of {L^). These coefficients are shown in 
Figure [5] for four voltages. The nonzero mean value of 
the angular momentum is a consequence of the asymme- 
try of these distributions with respect to m = 0. For 
example, the last two graphs at voltages U — 0.75 V and 
U = 1.0 V have a small asymmetry (about 1 per cent, 
not possible to see in the figure). This asymmetry will 
be much more pronounced in Model 2. 

At the beginning of this section, we discussed the role 
of symmetry in the models with respect to the transfor- 
mation Lp — >■ —ip and U — > —U. The parameters respon- 
sible for the asymmetry of the models with respect to 
these two transformations are the charged bridge poten- 
tial shift ipi and the right coupling shift Lpr- In Figure 
ini we show {Lz) as a function of fi and pr for Model 
Ic with the voltage fixed at [/ = 1.5 V; the red line in 
the figure marks the values actually used in Model Ic. 
The angular momentum first grows with pi but higher 
shifts suppress the effect again. This effect is similar to 
the Frank-Condon blockade for the current observed in 
Koch and von Oppen— . Selecting the optimum value for 
this parameter can enhance the angular momentum by a 
factor of three. The dependence of (L^) on ipr is shown 
in the lower graph of Figure |9l It is possible to maximize 
the " motor effect" by optimizing pr , gaining another fac- 
tor of 2. Coupling becomes symmetrical when (p^ = t^ ot: 
tpr — 27T as in Model lb, where (L^) vanishes. 



C. Results for the model 2 

We now discuss the more realistic model 2. Here 
the difference between the shapes of the potentials for 
the charged and uncharged molecule is large, leading to 
strong coupling and the impossibility to sort the vibra- 
tional transitions into more or less sharp lines. Further- 
more, we took a more realistic value for the moment 
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FIG. 9. Dependence of the mean value of angular momentum 
Lz at voltage U = 1.5 V on the model parameters of potential 
shift tfii (top) and coupling asymmetry i^a,. (bottom) in Model 
Ic. 



of inertia (corresponding to benzene), leading to a very 
dense vibrational spectrum. The calculation is therefore 
much more numerically demanding and we had to cover 
a smaller voltage range. 

The current and the excitation function of the molecule 
for three different temperatures are plotted in Figure [TU] 
as a function of voltage. As we have already mentioned, 
the density of the possible tunneling channels in Model 2 
is much higher than in Model 1. We therefore expect no 
distinct steps in the voltage dependencies of observables 
for this model. Another difference is that the ground 
state energies of the unoccupied and occupied bridge al- 
most coincide in Model 2. This means that the bridge is 
open for electrons with zero energy and tunneling events 
can happen at zero voltage. Despite the large amount of 
possible channels, we observe a zero current plateau in 
Fig.[TUl This is again a consequence of the small values of 
the Frank-Condon overlaps for the low-lying vibrational 
states in both potentials. The current plateau disappears 
for the calculation at higher temperatures, where higher 
states, not subjected to the Frank-Condon blockade, are 
already excited at zero voltage. 

Temperature effects become much more important for 
Model 2 than it was for Model 1. This is connected with 
the fact that even at the lowest temperatures considered 
here (50K), kT is comparable with distances between en- 
ergy levels of the molecule. It is also interesting to note, 
that the vibrational excitation curve in the lower graph of 
figure [TU] does not have the minimum at zero voltage, but 
at approximately U = ±0.15 V. This implies the pres- 
ence of a cooling effect of the current, which has already 
been observed experimentally^ and also been discussed 
theoretically2Si^i2^. 

Both the excitation function and the angular momen- 
tum voltage dependence (see Figure [TT]) are strongly 
asymmetric, which is not surprising, for the strongly 
asymmetric coupling. It is also obvious from Figure [3] 
that the left lead coupling strength is minimal in the 
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FIG. 10. The current and the excitation function for Model 
2a are plotted as a function of the bias voltage. The results 
for three different temperatures (including room temperature) 
are shown. The inset in the bottom part shows the contribu- 
tion to vibrational energy from the unoccupied and occupied 
bridge for the temperature 50 K. 



area where the unoccupied bridge wave functions arc lo- 
calized, which makes the overall coupling of the left lead 
smaller compared to the right lead. 

The angular momentum reaches much higher values 
(up to 12ft) inside the considered voltage interval, which 
could partially be attributed to the high moment of in- 
ertia. Temperature has some influence on the shape of 
the angular momentum curve but its maximum value is 
rather insensitive to temperature. From this we can con- 
clude that the motor effect can be observed at room tem- 
peratures as well as at cryogenic temperatures. 

As we discussed before, the nonzero mean values of 
the angular momentum are connected with asymmetries 
in the population distributions of the eigenstates of L^ , 
which are plotted in Figure [TU Large values of (L^) 
are accompanied by large asymmetries, which are clearly 
seen for U = -0.3 V and U = -0.4 V. 



V. CONCLUSIONS AND FUTURE PROSPECTS 

In this paper we have discussed the models of a 
molecule coupled to two conducting leads with coupling 
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FIG. 11. Mean value of the angular momentum of the 
molecule {Lz) as a function of voltage for Models 2a and 2b. 
The temperature dependence is also shown for Model 2a. 
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FIG. 12. Populations distributions pZi and p2 of the RDM 
in momentum basis are plotted for different bias voltages. It 
is clearly seen how distributions become wider and lose their 
symmetry while absolute values of the voltage grow. Graphs 
are plotted for room temperature (295 K). 



depending on the vibrational coordinate. It was shown 



how the steps in the current-voltage characteristics can 
be analyzed in the case of anharmonic vibrations of the 
molecule. We have also calculated the population of 
molecular vibrational states for different voltages. In the 
case of the model which was asymmetric with respect to 
an exchange of the left and the right leads, we had to dis- 
tinguish between the two charging states of the molecule, 
which have different populations, i. e. the population of 
vibrational states was strongly correlated with the charg- 
ing state of the molecule. 

In addition, we have studied the "motor effect", i. e. 
the response of the angular momentum of the molecule to 
the voltage applied across the junction. We have demon- 
strated that the mean value of the angular momentum 
strongly depends on the voltage. For the asymmetrically 
coupled molecule, the direction of the molecular rotations 
can be controlled by the polarity of the voltage. Signif- 
icant values of angular momentum were reached when 
the vibrational levels above the potential energy barrier 
against full rotation were populated. 

It will be interesting to see the effect of the higher or- 
der terms (see, for example, Esposito and Galperin— ), 
as well as friction, to the dynamics of the molecular rotor 
in future work. We believe that the "motor effect" must 
survive these additional corrections since it is quite stable 
with respect to changes to the model and it is a conse- 
quence of the breaking of symmetry of the junction with 
respect to the inversion of the angular coordinate and 
the interchange of the left and right electrode (chirality 
of the junction). 
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